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APPROXIMATING EXPLICITLY THE MEAN REVERTING CEV PROCESS 


N. HALIDIAS AND I. S. STAMATIOU 


Abstract. In this paper we want to exploit further the semi-discrete method appeared in Halidias and 
Stamatiou (2015). We are interested in the numerical solution of mean reverting CEV processes that appear 
in financial mathematics models and are described as non negative solutions of certain stochastic differential 
equations with sub-linear diffusion coefficients of the form ( xt ) q , where ^ < q < 1. Our goal is to construct 
explicit numerical schemes that preserve positivity. We prove convergence of the proposed SD scheme with 
rate depending on the parameter q. Furthermore, we verify our findings through numerical experiments and 
compare with other positivity preserving schemes. Finally, we show how to treat the whole two-dimensional 
stochastic volatility model, with instantaneous variance process given by the above mean reverting CEV 
process. 


1. Introduction. 


Consider the following stochastic models 


(i.i) 


St = So + fo M ■ S u du + / 0 4 (K) P • S u dW u , 

V t = V 0 + />! - k 2 V s )ds + f* k 3 (V s )idW s 


t&[0,T], 

t€[0,T], 


where St represents the underlying financially observable variable, V) is the instantaneous volatility when 
p = 1 or the instantaneous variance when p = 1/2 and the Wiener processes Wt, Wt have correlation p. 

We assume that V t is a mean reverting CEV process of the above form, with the coefficients ki > 0 for 
* = 1,2,3 and q > 1/2, since the process V) has to be non-negative. To be more precise the above restriction 
on q implies that V) is positive, i.e. 0 is unattainable, as well as non explosive, i.e. oo is unattainable, as can 
be verified by the Feller’s classification of boundaries 113 Prop. 5.22]. The steady-state level of V t is k\/k 2 
and the rate of mean reversion is k 2 . 

The system O) for p = q = 1/2 is the Heston model. When q = 1 we get the Brennan-Schwartz model 
pf] Sec. II], which apparent its simple form, cannot provide analytical expressions for St- 

Process V t for q = 1/2, also know as the CIR process [6] Rel 13], by the initials of the authors that 
proposed it for the term structure of interest rates, has received a lot of attention and we just mention two 
latest contributions to the study of such processes (see [li, jT2] and references therein). 

Process Vfforl/2<g<l has been also considered for the dynamics of the short-term interest rate [5j 
Rel (1)]. The stationary distribution of the process has also been derived in [2] Prop 2.2]. 

We aim for a positive preserving scheme for the process Vt. The scheme which we propose, and denote 
it semi-discrete (SD), preserves this analytical property of V) staying positive. The explicit Euler scheme 
fails to preserve positivity, as well as the standard Milstein scheme. We intend to apply the semi-discrete 
method for the numerical approximation of Vt in model (ED with 1/2 < q < 1 and compare with other 
positivity preserving methods such as the balanced implicit method (BIM) (introduced by [23 Rel (3.2)1 
with the positive preserving property fl6l Sec. 5]) and the balanced Milstein method (BMM) [16] Th. 5.9]LH 
Finally, we approximate the stochastic volatility model of (11,11) with p = 1/2. In jT5] a thorough treatment 
can be found, where also another stochastic volatility model is suggested. 

Section [2] provides the setting and the main results, Thcorems l2.ll and l2.21 concerning the C 2 — convergence 
of the proposed Semi Discrete (SD) method to the true solution of mean reverting CEV processes of the 
form of the stochastic volatility in (11.11) , as well as Theorem 12.31 which concerns the analogues of Theorems 
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^We give in the Appendix the form of all the above schemes for the approximation of Vt. 


1 



2 


N. HALIDIAS AND I. S. STAMATIOU 


m and 12.21 following an alternative approach. The main ingredient of this approach, inspired by m , is the 
simplification of the numerical scheme proposed, by altering the initial Brownian motion ( Wt ) to another 
Brownian motion (Wt) justified by Levy’s martingale characterization of Brownian motion, yielding to the 
same logarithmic rate as in Theorem 12.II but to a better polynomial rate \(q — instead of \(q — A | 
as shown in Theorem 12.21 

Section [3] is devoted to the proof of Theorem 12.11 while Section 0] and [5] concern the proofs of Theorems 
12.21 and l?~dl respectively. Finally, Section [G] presents illustrative figures where the behavior of the proposed 
scheme, regarding the order of convergence, is shown and a comparison with BIM and BMM schemes is 
given. In Section [7] we treat the full model (11.11) for a special case. Concluding remarks are in Section [8] 
and in Appendix [A] we briefly present numerical schemes for the integration of the variance-volatility process 

(Vt). 


2. The setting and the main results. 

We consider the following SDE 

(2.1) x t =x 0 + f (ki-k 2 x s )ds+ f k 3 (x s ) q dW s , te[0,T], 

Jo Jo 

where k±, k 2 ,k 3 are positive and 1/2 < q < 1. Then, Feller’s test implies that there is a unique strong solution 
such that xt > 0 a.s. when xo > 0 a.s. Let 


( 2 . 2 ) 


fe(x,y) = ki - k 2 ( 1 - 9)x - 




4(1 + k 2 e A) 


v 2q ~ l — k 2 0y + 


ki 


4(1 + k 2 6 A)' 

h(x) 


r 2 <?—i 


h(x,y) 

and 

( 2 . 3 ) g(x,y) = k 3 x q ~?^/y, 

where f(x,x) = a(x) = ki — k 2 x and g(x,x) = b(x ) = k 3 x q . 

Let the partition 0 = to < t± < ... < tw = T with A = T/N and consider the following process 

yf D (q) = yt n + h(yt n ,y t ) ■ a + f 2 (yt n )ds + J^ sgn(z s )g(y tn ,y s )dW a , 

with i/o = xq a.s. or more explicitly 


yt D (d) = Vt n + (ki- k 2 ( 1 - 0)y t „ - 


ki 


4(1 + k 2 e A) 


(yt n ) 


2q ~ 1 - k 2 9y^j ■ A 


ki 


-(ytj 29 1 ds + k 3 (y tn ) q 2 sgn(z s )y/yIdW s 


(2 ' 4) Jt n 4(1 + k 2 9A) Jt 

for t € (t n ,t n+ 1 ], where 9 € [0,1] represents the level of implicitness and 

k 3 


(2.5) 
with 

( 2 . 6 ) 




2 In ■= ytn 


1- 


k 9 A 


k\ A 


ki 


A- 


l + k 2 9Aj ‘ l + k 2 9A 4(1 + k 2 9A) 2 

Process (IQl is well defined when y n > 0 and this is true when ( 1+fc 1 2 g A ) fcl < 4(fc 2 Afci) and A(2 — 9) < 
Furthermore, (12.41) has jumps at nodes t n . Solving for y t , we end up with the following explicit scheme 


(2.7) yf D (q)=Vn+J 


k 2 


r (yt n ) 2q ~ 1 ds + 


k 3 


4(1 + k 2 9A) 2 ' 9 ‘ n/ ““ ' 1 + k 2 9A 

with solution in each step given by m Rel- (4-39), p. 123] 

y? D ( q ) = (z t ) 2 , 

which has the pleasant feature yf D (q) > 0. 


(ytn) 9 


sgn(z s )y/yldW s 
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Assumption A Let the parameters ki,k 2 ,k 3 be positive and such that ( 1+fc 1 2 g^) fcf < 4(fc 2 A fci) and 
consider A > 0 such that A(2 — 6) < for 9 £ [0, 1 ]. Moreover assume xq > 0 a.s. and E(xo) p < A 
for some p > 4. 

Theorem 2.1. [Logarithmic rate of convergence] Let Assumption A hold. The semi-discrete scheme wa 
converges to the true solution of \2.1\) in the mean square sense with rate given by 

( 2 -8) E sup \yf D {q) - x t \ 2 < —===, 

o <t<T 

where C is independent of A and particularly 

C := 72y^(fc 3 ) 4 T 2 e 6T2fe2+fc2T , 

where e is such that 0 < e < q A (q — j). 

Assumption B Let Assumption A hold where now xq £ R and Xq > 0. 

Theorem 2.2. [Polynomial rate of convergence] Let Assumption B hold. The semi-discrete scheme \2. 7[ ) 
converges to the true solution of \2.1\) in the mean square sense with rate given by, 

(2.9) E sup I y? D (q) - x t \ 2 < CA^~^ A i, 

0 <t<T 

where 

C := 8 ^12fc|T(^A 4g (xo + hT^CkM.A V 1) \J 3k 2 3 T^/A 2 (x 0 + fciT) 2 f A 4q _ 2 
x 

and Chk is a constant described in j4-ll\ ) and A is an appropriately chosen positive parameter which satisfies 
\4-Bty and always exist, p(A) := 2 (i - q yi(k^ yz — 1> quantity Ck 2 ,k 3 , 0 ,a is given in Lemma \3A\ and e > 1. 



Inspired by El we remove the term sgn(z s ) from (12.41) by considering the process 


( 2 . 10 ) 


W t = sgn(z s )dW s , 


which is a martingale with quadratic variation < Wt, Wt >= t and thus a standard Brownian motion w.r.t. 
its own filtration, justified by Levy’s theorem [T7l Th. 3.16, p.157]. Therefore, the compact form of (12.41) 
becomes 

Vt D = x o + [ ( fc i - & 2 (1 “ - k 2 Oys) ds 

Jo 

+ jt (ki- k 2 { 1 - 9)y tn - nf q ~ X ~ ki9yfj ds + k 3 s{y s ) q ~^ ^y[dW s , 

for t £ (t„,t n +i]. Consider also the process 

(2.11) x t = x 0 + / (fci - k 2 x s )ds + 

Jo 

The process (x t ) of (12.11) and the process (x t ) of (12.111) have the same distribution. We show in the following 
that Esup 0<t<r | yt D (q) — x t | 2 —t 0 as A J. 0 thus the same holds for the unique solution of (12.11) . i.e. 
Esup 0<t<T | yf D (q) — Xt\ 2 —> 0 as A j. 0. To simplify notation we write W, (x t ) as W, (x t ). We end up with 
the following explicit scheme 


k 3 (x s ) q dW s , t £ [0, T]. 


( 2 . 12 ) 


v? (?) = yn + 


k'i 


4(1 + k 2 9A)‘ 


: (Vt n ) 2q ~ 1 ds + 


1 + k 2 9A 


(y tn ) q 


y/[hdW s 


where y n is as in (EH). 
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Theorem 2.3. [Logarithmic and Polynomial rate of convergence] 

Let Assumption A hold. The semi-discrete scheme \2.12\ ) converges to the true solution of \2.1 )) in the 
mean square sense with rate given by 

(2-13) E sup \yf D (q) — x t \ 2 < —====, 

o <t<T ym(A) _i 

where C is independent of A and given by 

C := 32y^(fc 3 ) 4 T 2 e 6r2fe2+fc2T , 

where e is such that 0 < e < q — ~. 

In case Assumption B holds, the semi-discrete scheme \2.12\) converges to the true solution of \2.1\) in the 
mean square sense with rate given by, 

(2.14) E sup \yf D (q) -x t \ 2 < CA^ q ~i\ 

0 <t<T 

where 

C := y (l2klT(y/A 4q (x o + k 1 T)^C k2tk3! g A V 1) \/ 3k 2 3 T^A 2 (xo + A Aq _fj 

x ^e 6fe2T2 + ^|(xo) (1 - 9) " (A) 

and Chk is the constant described in and A is an appropriately chosen positive parameter which 

satisfies J 4-12\ ) and always exist, u(\) := 2 (i-<p i (k 3 )‘ i ~ 1) Quantity Ck 2 ,k 3 ,g ,A is given in Lemma \3-4\ and 
e > 1. 

In the following sections we write for simplicity y^ D or y t for yf D (q). 


3. Logarithmic rate of convergence. 
We rewrite (12.41) in a compact form 


Vt D = xo + [ (ki- k 2 ( 1 - 0)y s - k 2 dyg) ds 

Jo 


+ [ [ fci - fc 2 (l ~ 0)y tn - , l flA' feJ 211 1 ~ k 2 dy t ) ds + fc 3 f sgn(z s )(y s ) q ^y/yfdW s , 


4(1 + k 2 e A) 

for t £ ( t n ,t n+ \ ] where 

s = t s£(t t -iil 7=0 n s=I tj+1 ’ for s e [tptj+il a _ o 
s t 3 ,se(t 3 ,l 3+ i|,j u ,...,n, 1 for s £ , t] J “ ’' ’' 


n — 1. 


3.1. Moment bounds. 

Lemma 3.1 (Moment bound for SD approximation). It holds that 

E sup (y t ) p < ApEixo + k^P, 

0 <t<T 

for any p > 2, where A p := exp j gdriD fcg ^p-j_ + 2L±^ t }. 

Proof of Lemma 1X71 We first observe that ( y t ) is bounded in the following way 

f>t /»t i 

0 < y t < x 0 + kids + k 3 / sgn(z s )(y s ) q ~^ 3 Jy]dW s 


JO 


JO 


< Xo + hT + ks sgn(z s )(y§) q = y/y[dW s := u t 


Jo 

a.s., where the lower bound comes from the construction of ( y t ) and the upper bound follows from a com¬ 
parison theorem. We will bound ( Ut ) and therefore (yt), since 0 < yt < Ut a.s. Set the stopping time 
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tr := inf{t £ [0,T] : u t > R}, for R > 0 with the convention inf 0 = oo. Application of Ito’s formula on 
i u tAr R Y implies 


(ut A T R ) p = {x 0 + k x T) p + p(j> 2 J ( u s ) p 2 (ys) 2q 1 y s ds 

rtAT R 

+pk 3 / sgn{z s ){u 8 ) p ^ 1 {ys) q ~^^/yldW s 
Jo 


viv-i) r tATR 

< ( Xo + k 1 T) p + PKP 2 jf (u s ) p ~ 1 (ys) 2q ~ 1 ds + M t 


< ( Xq + kiT) p + 


< {x 0 + k\T) p + 


Pip ~ 1) [ tATR ( P~ 1 

2 P 

PiP ~ !), 2 fP - 1 , 2 P_ 


2 P 


-1 


( Us )P +- {y s )^-^\ d s + Mi 


-l \ r tAT R 


2 P P 


(u s ) p ds + M t , 


where in the second step we have used that 0 < y t < u t , in the third step the inequality x p l y < e^y-x p + 
p_i y p , valid for x A y > 0 and p > 1 with e = |, in the final step the fact \ < q < 1 and M t := 

fg ATR sgn(z s )(u s ) p ~ 1 (y 8 ) q ~^^/yfdW s . Taking expectations in the above inequality and using that M t is 
a local martingale vanishing at 0, we get 


E(ut A .J p < E(x 0 + k 1 T) p + PiP 2 + 


< E(xo + k\T) p exp 

< ApE(xg + kiT) p , 


pip -0,2 fp- 1 , 2 P “ 


+ 


2 P P 


E(u sA r R ) P ds 


T 


where we have applied Gronwall inequality jH Rel. 7]. We have that 

iVtA t r T = iyr R ) P l { r R <t } + iytn { t<r R} > iyt) P I { t<r R }, 

thus taking expectations in the above inequality and using the estimated upper bound for E(it tATil ) p we 
arrive at 

E iyt) p l{t<r R } < E(ytAr R ) p < E (u tATR ) p < A p E(x 0 + k^Y, 
and taking the limit as R —> oo, we get 

lim E(y t Yl{ t <r R } < A p E(x 0 + fciT) p . 

R—too 

Let us fix t. The sequence of stopping times tr is increasing in R and t A tr — > t as R —> oo, thus 
the sequence iytYht<T R } is nondecreasing in R and {yt) p I{t< Tfl } —> ( y t Y as R —> oo. Application of the 
monotone convergence theorem implies 


(3.1) 


HytT < a p e(x 0 + hTY, 


for any p > 2. Using again Ito’s formula on (u t ) p , taking the supremum and then using Doob’s martingale 
inequality on the diffusion term we bound Esup 0<t<T (it* ) p and thus Esup 0<t<r (?/t) p . □ 

Lemma 3.2 (Error bound for SD scheme). Let n s integer such that s € [tn s j£n s +i]- It holds that 

E|j/ S - y~s\ P < 0 p A pl2 , E|y s - y 8 \ p < A p A p / 2 , 
for any p > 0, where the positive quantities A p ,A p do not depend on A. 
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Proof of Lemma kOl First we take a p > 2. It holds that 


\y s -ys\ p = / (fci - fc 2 (l - 0)yu - k 2 dyu)du + 


k20ygdu — 


k20y s du 


■n 


+j (^ - fc 2 (i - - 

5 P_1 




-(ytn s ) 2q 1 ^jdu + k 3 J^ sgn{z u )(yu) q ^y/y^dW^ 


4(1 + k 2 0A) 

(j J (fci - fc 2 (l - 0 ) 2 /* - fc 2 6M) du| P + fc 2 0 p (M P (M+i - ts;) p + fc 2 0 P M) P (M+i - s) p 


fci - fc 2 (l - % tns - 


-(^J 29 ” 1 du 


4(1 + fc 2 6»A) 

< 5P- 1 (its; - sr 1 jf |fc! - Ml - %a - ^0y s | p dw + fc p 0 p ((p s ~) p + (y s )P) A p 


fcf| f sgn{z u ){yu) q *y/fh.dW u \fj 


h - fc 2 (l - 6')y tn - 


fci 


4(1 + k 2 0A) 


(yt n .) 


2 <?-l 


A p -(- fc p | J" sgn(z u )(ya) q *y/y^dW u \ p f 


where we have used Cauchy-Schwarz inequality. Taking expectations in the above inequality and using 
Lemma O and Doob’s martingale inequality on the diffusion term we conclude 

(3.2) E|y s - y s | p < i p A p / 2 , 

where the positive quantity A p except on p, depends also on the parameters M fc 2 , M 9, q, but not on A. 
Now, for 0 < p < 2 we get 

% s - ys \ p < (%s - ys \ 2 ) p/2 < i P A p/2 , 

where we have used Jensen inequality for the concave function 4>{x) = x p ' 2 . Following the same lines, we 
can show that 

(3.3) E|y s - < I P A p/2 , 

for any 0 < p, where the positive quantity A p except on p, depends also on the parameters M M fc 3 , 9, q, 
but not on A. 

For the rest of this section we rewrite again the compact form of (12.41) in the following way 


□ 


(3.4) 


Vt = x o+ / fe{ys,ys)ds+ / sgn(z s )g(y s ,y s )dW s + / fi{yt n ,yt)ds, 


where is given by (12.21) and the auxiliary process (h t ) is close to (y t ) as shown in the next result. 

Lemma 3.3 (Moment bounds involving the auxiliary process). For any s £ [0, T] it holds that 

(3.5) E| h s - y s | p < C p A p , E\h s \ p < C h , 
and for s £ [t n ,t n + 1 ] we have that 

(3.6) E\h s -y s \ p <CpA p ' 2 , E\h s — ys\ p < C p A p k 2 , 
for any p > 0, where the positive quantities U p , Cp> C p , Ch do not depend on A. 

Proof of Lemma \3.A We have that 


I h s -y s \ p = 


h(yt n ,yt)du 


< \t n+1 - s\ p \h(y tn ,yt)\ p , 


for any p > 0, where we have used (i3~H) . Using Lemma ET! we get the left part of (13.51) . Now for p > 2 and 
noting that 

E\h s \ p < 2 p - 1 E\h s -y s \ p + 2 p ~ 1 E\y s \ p 

< 2 P ~ 1 C P A P + 2 p -M p E(x 0 + k\T) p < C h , 
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we get the right part of (13.51) . where we have used Lemma l3Tl The case 0 < p < 2 follows by Jensen’s 
inequality as in Lemma 13.21 

Furthermore, for s £ [f n ,f n .|_i] and p > 2 it holds 

nK-y s \P < 2P- 1 E\h s -y s \P + 2P- 1 E\y s -y s \P 
< 2 P ~ 1 C P A P + 2P- l A p A p / 2 < C p A p / 2 
where we have used m and in the same manner 

E\h a - yj\ p < 2 P ~ 1 C P A P + 2 p - 1 A p A p/2 < C p A p / 2 . 

The case 0 < p < 2 follows by Jensen’s inequality. □ 


3.2. Convergence of the auxiliary process [h t ) to (xt) in C 1 . We first estimate the probability of z t 
being negative when at the same time y tn > A 1-2 ^, for 0 < £ < 4. 

Lemma 3.4. For every t £ [t n ,tn+i] it holds 

(3.7) P ({** < 0} n {y tn > A 1 - 2 ?}) < C kMA VA, 

where Ck 2 ,k 3 ,6 A '■= ^ g ^ A and A(2 — 8) < and ( - 1+ ^ 6 , A ) < 4 k 2 . Relation \3. 7| ) implies that 
p ({zt < 0} n {y tn > A 1 -^}) = 0(V A), as A f 0. 

Proof of Lemma \3.4\ By the definition (12.51) of (z t ) for t £ [t n ,tn+ 1 ] and for 0 < f < |, we have that 
A := {zt < 0} D {yt n > A 1 -^} = | (yt n ) q -Hw t - W u ) < - 2(1 \ k f A) Vy^} n{y tn > A 1 -^} 


(3.8) 

where 

and 


A i := < Wt — W, < - 


C Ai U A 2 , 

2(1 + k 2 6A) 


A 2 :={W t - W. < - 


Vy^(ytJ- q+h ]n{ ytn > 1 }, 

2(1 \ k f A]] Vy^(yt n )~ q+i } n {1 > y tn > a 1 ^}. 


The following inclusion relations hold for the event Ai, 
Ai 

2(1 + MA) 


C AW„ < — 


I ytn 1 - 


(k 3 ) 2 a 


k 2 A 

1 + k 2 9A J 4(1 + k 2 9Ay 


(ytJ 29_1 (2/tJ 9+ " } n{y tn > 1 } 


c 


c 


AW n < -- 
AW„ 


k 3 


/1 — fc 2 (2 — 9) A 

(k 3 ) 2 A } 

j 1 + k 2 9A 

4(1 + MA) 2 J 


< - 


2 ^/(l — k 2 (2 — 9)A)(1 + k 2 9A) 


\ft> tji k 3 \ft 1 tri 

1 ( fc 3) 2 


when A(2 — 9) < and ( 1 _ h fc 2 g A ) < 4 k 2 , where A W n := W t — W tn . It holds that 


(3.9) 


P(G < -p) = 


1 


1 


-.e~ u2 / 2 du < / e- u2 ' 2 du = / e~ u2 / 2 du < t e ~{f>Y/ 2 , 
J —oo J B P 


/—oo V / 27 t 

for every standard normal random variable G, where in the last step we have used [T7[ Ineq. (9.20), p.112] 
valid for /3 > 0. Using the fact that -^='= is a standard normal r.v. and ignoring the exponential term in 
(EH), since its exponent is negative, we get that 

k 3 


2^{\-k 2 {2~9)A) 


y/t tfi A Gfc 2 & V^A. 


(3.10) 


P(Ar) < 



























8 


N. HALIDIAS AND I. S. STAMATIOU 


The following inclusion relations hold for the event A 2 , 

A 2 


C { A W n < - 


2(1 + k 2 9A) / 1 — k 2 (l — 9)A faA 


fa 


iyt„ 


+ 


1 + k 2 9A 1 + faOA 4(1 




n{l > y tn > A 1 - 24 } 


C { AW n < - 


2(1 + k 2 9A) l A1 _ 2i l-k 2 (l-6)A 


fa 


1 + fad A + \ kl 4(1 + k 2 6A)) 1 + k 2 0A f 


( fa ) 2 


A } 


c 


AWn < _2_ ya^faiT^jAxrTfaMj aI _ £ 

\Jt t n k 2 yjt t n 

(fcs) 2 


when A(1 — 6) < W an d (1 _j_ fc 3 2 g A ^ < 4fci- Using again (13.91) we have that 

fa 


P(A 2 ) < 


(3.11) 


< 


2^/(1-A 2 (1-0)A) 
fa 


_ , , 2 II IJ 2 ll-g)A)(l + lj 2 eA) A l-2£ 

A«-5e UW 


A^e ( fe 3 ) : 


r(l-fc 2 (l-e)A)(l+fc 2 eA)A~ 2 « 


2^/(1-fe(l-0)A) 

Taking probabilities in the inclusion relation (13.81) and using (13.101) and (13.111) we get 
P(A) < P (A 1 ) + E(A 2 ) 


< Cfe 2 ,fc 3 ,e, aVA + 

< C'fe 2i fe 3 ^,AV / A, 


fa 


2\/ (1 — fc 2 (l — &)A) 


Aie -j^(l-k 2 (l-8)A)(l+k 2 eA)A- 2 S 


since A^e A 25 = o(VA) as A j. 0. Finally, note that Ck 2 ,k 3 ,e,A = — 

y'l-/c 2 (2-0)A 


the 0(\/A) notation, (see for example [24] 1. 

We will use the representation (13.41) and write 

(3.12) h t -x t = (fe(ys,ys) ~ fe(x 8 ,x s )) ds + / (sgn(z s )g(y s ,y s ) - g(x s ,x a )) dW s . 

Jo Jo 

Proposition 3.5. Let Assumption A hold. Then we have 


fa as A 4. 0 which justifies 
□ 


(3.13) 


, , / Aj A«-5 „ 1 \ , T 

sup E\h t - x t \ < J 2 -h J3-h 3fc|T— e k2 , 

0 <t<T \ rne m me m m 


for any m > 1, where e m = e m ( m + 1 )/ 2 an & 

J 2 := 12 klT(^JA<± q E{xQ + k\T) A( iCk 2 ,k 3 ,e,A V 1), J3 := 3k%T A 2 E(xq + k\T) 2 \JA^ q - 2 . 

Proof of Provosition \3.5\ Let the non increasing sequence {e m } m6 N with e m = e ~ m ( rn + 1 )/ 2 an d eQ = 1. We 
introduce the following sequence of smooth approximations of |ar|, (method of Yamada and Watanabe, [27] ) 

r M rv 

<t>m(x) = dy ip m (u)du, 

Jo Jo 

where the existence of the continuous function tj) m (u) with 0 < if m (u) < 2/(mu) and support in (e m ,e m _i) 
is justified by / e em_1 (du/u) = m. The following relations hold for <j> m £ C 2 (R,M) with <j> m ( 0) = 0, 

\x\ - e m _i < (f m (x) < \x\, \4>' m (x)\ < 1, iSK, 


We have that 
(3.14) 


1 4*rn (^0I < —rTi when e m < III < e m _i and \4>'L(x)\ = 0 otherwise. 
m\x\ 


E\h t - x t \ < e m -1 + E4> m (ht - x t ). 
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It holds that 


(3.15) 

and 


(3.16) 


fe(y$, Us) - fe(x a ,x a ) = (fci - fc 2 (l - 6)y a - k 2 0y a ) - (ki - k 2 x s ) 

= —fc 2 (l - 0)(ys - x 8 ) - k 2 0{y a - x a ) 

= k 2 {l - 9)(h s - y s ) + k 2 0(h s - yj) - k 2 (h s - x s ) 

\sgn{z s )g(y s ,y s ) - g(x s ,x s )\ 2 = \sgn{z s )k 3 (yg) q ~^^- k 3 (x s ) q \ 2 

< ^3 {iys) q ~^VTTs{sgn(z s ) -1) + Jyl ((yi) 9- ^ - (y s ) 9- ^) + ((y s ) q - (x s ) 9 )) 

< ^(y s -) 29_1 y s (syn(z s ) - l) 2 + y s ((y s -) 9- ^ - (y s ) 9- *) + (( y s ) q - (a; s ) 9 ) 2 ^ 

< 3fc 2 (( ys) 2q ~ 1 y s (sgn(z s ) - l) 2 + y s |y s - - y s | 29-1 + (J|y s - z s |) 2 ) 

< 3fc 2 ((yj) 29_1 y s (sgn(z s ) - l) 2 + y s |y s - - y a | 29-1 + | h s - y s \ + \h s - x s |) , 


where we have used properties of Holder continuous functions and namely the fact that x q is q —Holder 
continuous for q < 1, i.e. \x q — y q \ < \x — y\ q , and that x q is 1/2—Holder continuous since q > 1/2. 
Application of Ito’s formula to the sequence {</> m } m 6N, implies 


< 


(j)m(h t - x t ) = / fimihs - x s )(fe(ys,ya) - fo(x s ,x s ))ds + M t 
Jo 

i j 

+ 2 J - x s ){sgn{z s )g{y ai ys) - g{x s ,x s )) 2 ds 

[ (k 2 ( 1 - 6)\h s - y a | + k 2 0\h s - y a | + k 2 \h s - x«|) ds + M t 
Jo 


+3^3 


i|/l s x s 


(( ys) 2q 1 y s \sgn(z s )-l\ 2 + y s \ys-y s \ 2q 1 + \h s - y s \ + \h s - x s |) ds 


f f 3 k z f f 

< k 2 (l — 0) \h a - ys\ds + k 2 d / \h s - y a \ds-\ - — \h s - y s \ds + k 2 / |/i s - x s |ds + M t 

J o do me m y 0 do 

+ _^3_ f (y s ) 2 q- 1 y s \ Sgn ( Zs ) _ l| 2 rf S + -^3- f y s |y S - y s | 29_1 rfs + 

me m do me m do ™ 

where in the second step we have used (13.1511 and (13.161) and the properties of (j> m and 

M t := / <t>' m (h u - x u )(sgn(z u )g(yu,yu) - g(x u ,x u ))dW u . 
do 

Taking expectations in the above inequality yields 

r f rt 3k 2 r t 

E(j) m (h t - x t ) < k 2 (l - 9) E\h s -ys\ds + k 2 9 E\h s - y 8 \ds -\ - — E\h a - y a \ds 

do do me m do 

ft 51.2 ft 3k 2 T /■* 

■ / E(y s -) 29_1 y s |syn(2; s ) - l| 2 dsH-— / Ey s |y s - - y s | 29_1 ds + —— + k 2 / E|/i s -z s |ds 

m Jo mCm Jo m Jo 


3k 2 ^ 


< 


**(1 - 0)TCi^A + MTCiVA + -®——A + /c 2 / E|/i a -x a |d. 

TTlGm J 0 


0 

3fcs ^ E(i/s) 29_1 y s |syn(2; s ) - l| 2 ds + [ jE(y s ) 2 \/E\ys - y s | 4 «“ 2 ds + 


men 


jo 

< fc 2 r((i - 0)61 + 0 Ci)Va + 


me. 


m jo 


m 


3k 2 TC\ 

me m 


A + k 2 E|h s — x s |ds 


3fc: 


2 /•* 


E( yg) 2q 1 y s \sgn(z s ) 1| z ds + 


me. 


m JO 


3kiT 


me r 


2 A/i 4 o _ 2 A 9 -^ + ®, 


m 
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where we have used Lemma l3.3[in the second step and Holder inequality, Lemmata 13.II and 13.21 in the third 
step and the fact that E M t = 0o It holds for s €E [t n ,t n + 1 ] that 

E ( ys) 2q ~ 1 y s \sgn(z s ) — 11 2 = E (4(j/ t J 2?_1 2/ s I{z.<°}) 

< 4E| {y tn f q ~ l ys - {yt n f q \ + 4E ((y t + 4E ((Vt n ) 2q I { , m < 0} I{y tn >Ai-* } ) 

< 4E \ (jjt n ) 2q ~ 1 {ys ~yt n ) | + 4A 2 «- 4 «f + 4y/E(y t J 4 «y/F ({* s < 0 } n {y tn > A 1 ’ 2 ?}) 

< 4^E(y t J 4 ?-V E l2/ s - 2/d 2 + 4A 2 «" 4 ^ + 4^^ C kMA yfZ 

< ^A 4q _ 2 E(x 0 + + 4A 2 «- 4 ^ + 4 v / H 4g E(x 0 + k^y/C^^Ai, 

where we have used Holder inequality and Lemma 13.41 in the third step and Lemmata 13.21 and 13.II in the final 
step. For £ = \ — yfo we get the estimate 

(3.17) E(y s ) 2q ~ 1 y s \sgn(z s ) - 1| 2 < 4 ^A iq E(x 0 + k ± T) 4 iC k2 ,k 3 ,e,A V 1 ^ A*. 

Thus (13.141) becomes 

E\h t — x t \ < e m _i + J\\fK + 3kiTC\ - V Ji -b J 3 -b 3fc|T- k 2 [ E|/i s — zJrfs 

me m rne m me m m ,/ n 


< J 2 


A q 


1 


men 


J 3 - 1 -3fcoT-b k 2 / E|/i s — * s |ds 

me m to J 0 


< J2- 


J 3 


A q ~ 2 


3fc|T— 
mem to 


3^2 t 


where in the second step we have used the asymptotic relations, A K = o( A 9 ~5) as A | 0 for any k > 
1/2, e m _i = o(A') as to -A 00, vA = o(^—) for any re < 1 as to —00, in the last step we have used the 
Gronwall inequality and J 3 is as defined in Proposition 13.51 while 

Ji ■- k 2 T({i-e)c 1 + ec 1 ). 

Taking the supremum over all 0 < t < T gives (13.131) . □ 

3.3. Convergence of the auxiliary process (h t ) to ( xt ) in C 2 . 

Proposition 3.6. Let Assumption A hold. Then we have 

E sup \h t — x t \ 2 < , 

0 <t<T yln(A ) -1 


(3.18) 

where C e is independent of A and given by C e := 72 A(ko)^T 2 e 6T fc2 + fc2T , where e is such that 0 < e < 
I A (5- !)• 

Proof of Provosition VTT\ We estimate the difference \£ t \ 2 := \h t — x t \ 2 . It holds that 

\£t\ 2 = / if{y§,yz) ~ f(x s ,x s )) ds + / (sgn(z s )g(y s ,y s ) - g(x s ,x s )) dW s 

Jo Jo 

<2 T f (k 2 (l - 9)\h s - y s \ + k 2 d\h s - ys\ + k 2 \S s \) 2 ds + 2\M t \ 2 

Jo 

< 6T/c|(l — 0) 2 f \h s -y s \ 2 ds + 6Tk 2 d 2 [ \h s - y 7 \ 2 ds + 6Tk 2 [ \S s \ 2 ds + 2\M t \ 2 , 

J 0 Jo Jo 

where in the second step we have used Cauchy-Schwarz inequality and (|3.15p and 


M t := / {sgn{z u )g{y^yu) - g(x u ,x u ))dW u . 
Jo 


2 The function d(u) = (^^(hu — x u )sgn(z u ){g(yu, Vu) — g{x u , x u )) belongs to the space -M 2 ([0, £];R) of real valued measurable 

Tt~ adapted processes such that E fg \d(u)\ 2 du < oo thus 1211 , Th. 1.5.8] implies E Mt = 0. 
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Taking the supremum over all t G [0,T] and then expectations we have 

E sup \£ t \ 2 < 6Tfcf(l - Of [ E|/i a -y a \ 2 ds + 6Tk%0 2 [ E\h s -y 7 \ 2 ds 
0 <t<T Jo Jo 

+6 Tk 2 I E sup \£i\ 2 ds+ 2E sup \M t \ 2 
Jo C KKs o <t<T 


< 6 T 2 mi-eyA 2 A + 6 T 2 k^ 8 2 A 2 A + 6 Tk 2 / E sup |£;| 2 ds + 8E|M T | 2 , 

Jo 0<l<s 


(3.19) 


where in the second step we have used Lemma l3.2l and Doob’s martingale inequality with p = 2, since M t is 
an R—valued martingale that belongs to C 2 . It holds that 


E|M t | 2 := E 


\sgn(z s )g(y$,y s ) - g(x s ,x s )\dW s 


= E 


\sgn(z a )g(y s ,y a ) - g(x s , x s )\ 2 ds 


< 3/c|E^y (( y s ) 2q 1 y s (sgn{z s ) - l) 2 + y s \ys - y s \ 2q 1 + \h a - y s \ + \h a - x g |) ds^j 

< 3k\ f E ((ys) 2q ~ 1 y s \sgn(z s ) - 1| 2 ) ds + 3k 2 f E (y a \y$ — y s | 29-1 ) ds 

Jo Jo 

+3 k 3 I E\h s — y a \ds + 3k 2 I E|£ s |ds, 

Jo Jo 

where we have used (13.1611 . Now, we use again the estimate (13.171) to get 

E\M t \ 2 < J 2 AJ + J 6 VA 2 ?- 1 + 3k 2 TCiA + 3 kj [ E\£ a \ds 

Jo 

< J 2 Ai + J 6 A 9 -5 +3k 2 [ E\£ a \ds, 

Jo 

where we have used the asymptotic relations, A* = o(A 9- ^) for all l > i as A 4- 0, the quantity Jq is given 

by J& := 3k 3 Ty / A 2 E(xo + k\T) 2 ^j ^ 45-2 and J 2 is as given in the statement of Proposition [33] and depends 
on A through Ck 2 ,k 3 ,e,A where as already stated before we have that Ck 2 ,k 3 ,e ,A —> k 3 , as A 4 . 0. 

Relation (13.191) becomes 

rT rT 

E sup |£ t | 2 < 8 J 2 A 3 + 8 JeA q ~ 2 + J 5 A + 6Tfc 2 / E sup |£z| 2 ds + 24fc 2 / E sup \£i\ds 
0 <t<T Jo 0 <l<s Jo 0<Ks 


A 


A 9- 


1 


< 8J 2 A3 + 8 J 6 A ci ~i + 2AkiT J 2 -+ J 3 -+ 3 k 2 T— e k2± + 6 Tk 2 / E sup \£i\ 2 ds 

Y me m me m mJ J 0 0 <i<s 

< 24/c?TJ 2 e fe2T+6T2fe2 ^— + 24klTJ 3 e k2T+6T2k2 ^— + 72 (fc 3 ) 4 T 2 e fe2T+6T2fe2 —, 

me m rne m m 

where we have used Proposition 13.51 in the second step with the sequence e m as defined there, Gronwall 
inequality in the last step and the asymptotic relation A K = o( £ e ) as m —> 00 , for any k > 0 and J 5 is 
independent of A and given by J 5 := 6 T 2 fc 2 ((l — 0) 2 A 2 + 0 2 A 2 ). 

We take m = Vln A _A , with A > 0 to be specified soon and note that e ' /lnA x = o(A^ A ), as A 4 . 0, since 
gvdn n _ 0 ( n ^ as n —> 00 . Moreover it holds that 


A «- 


A 9 "5 to A 9 "* ~ _1_3A e ^ VlnA - x 

-=- e 2 vlnA = A 9 ^ -2- 

m z _ In A~ A 

e 2 e 2 


A- 
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Now, since q > h there is an e > 0 small enough such that q — i — e > 0. We take A = ^ and conclude that 


A 9 2 i ,e 2 

-= A 9- 2 -e 


i\/lnA~TT 


~t 0, 


Cm A 3 

as A —> 0 which in turn implies the asymptotic relation ^ T = o(—) as m —> oo, with the logarithmic rate 

1 

stated before. In the same way we can show = o(^) asm-) oo, by taking e < We finally arrive at 

1 


E sup \£ t \ 2 <72{k 3 ) 4 T 2 e k2T+6T k2 - _, 

o <t<T VlnA"w 


by taking 0<e<^A(g - A), which implies (13.181) . 


□ 


3.4. Proof of Theorem 12.11 In order to finish the proof of Theorem 12.II we just use the triangle inequality, 
Lemma [3731 and Proposition 13.61 to get 


E sup \y t -x t \ < 2E sup \h t -y t \ +2E sup \£ t \ 

0<t<T 0 <t<T 0 <t<T 

a 

2 Z- ■ ■■ 


< 2C 2 A 2 + 2- 


< 


C 


Vln A -1 -\/ln A -1 ’ 

where C = C{k 2 , ko, e, T), and given in the statement of Theorem 12.II 

4. Polynomial rate of convergence. 

We work with the stochastic time change inspired by [3]. We define the process 

192(fc 3 )V 


7(f) := f 
Jo 


rds 


o [(y s y- q + (x s y- q Y 

and the stopping time 

r; := inf{s £ [0, T] : 6 Tk 2 s + 7 (s) > l}. 

The process 7 (f) is well defined since x t > 0 a.s. and y t > 0 (see Section 0). 

The difference \£ t \ 2 := \h t — x t \ 2 is estimated as in Section [3] and we get, as in (13.191) . that 


(4.1) 


E sup \£t \ 2 < J 5 A + 6 Tk 2 ( E sup \£i\ 2 ds + 8 E\M T \ 
0 <t<T Jo 0 <l<s 


where r a stopping time and J 5 independent of A is as in proof of Proposition 13.61 The main difference here 
will be the estimation of the last term in (HU)- The approach in Section [3] resulted in the C 1 estimation 
E|£t| where we used the Yamada-Watanabe approach. Now, we use the Berkaoui approach. Relation (13.161) 
becomes 

\sgn(z s )g{y s ,y s ) - g{x s ,x s )\ 2 < 3 k 2 (( ys) 2 q ~ 1 y s (sgn{z s ) - l ) 2 +y s \ys - 2/ s | 29-1 + | (y s ) q - (x s ) q | 2 ) 

< 3 k 2 ({ysy^yaisgnizs) - l ) 2 + y s \y s - j/ s | 29_1 ) + \(y s ) q - ( x s ) q \ 2 ((y s ) 1-9 + (a : s ) 1_9 ) 2 


< 3 kl((y s ) 2q 1 y 3 (sgn(z s )-l) 2 + y s \y s -y s \ 2q x ) + g (\h s - y s \ 2 + \£ s \ 2 ) (is)', 

where we have used the inequality 

(4.2) |a 9 - 6 9 |(a 1-9 + 6 1 " 9 ) < 2q\a - b\, 

valid for all a > 0, b > 0 and \ < q< 1. Consequently, we get the upper bound 

pT 2 

E|M r | 2 :=E / \sgn(z s )g(ys,y s ) ~ g(x s ,x s )\dW s 
Jo 

< 8J 2 A* + J 6 A q ~i + ^ [ E\h s -y s \ 2 ( 7 s )'ds+^- [ E\S s \ 2 (^ s )'ds 

° Jo ° Jo 

< 8J 2 Ai + J 6 A q ~? + l [ y /E\hs-ys\WH('YaY) 2 d3+l / E| £ s \ 2 { ls )’ds, 

° Jo ° Jo 
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where we used the estimate (13.171) and Holder inequality J 2 is as in the statement of Proposition [331 and Jo 
independent of A is as in proof of Proposition 13.61 Relation (Id.II) becomes 

E sup \£ t \ 2 < 8J 2 A* + 8 J 6 A q ~^ + 6Tfc 2 [ E sup \£i\ 2 ds + [ y/M\h s - y s \*y/E{('y B )') 2 ds 
0 <t<r Jo 0 <l<s Jo 

+ [ E|£ s | 2 ( 7s )'ds 


< 8(J 2 VJ e )A(«-5) A 4^A 2 / * e( - 192(fc3)2<?2 - J] ds + f E sup \£i\ 2 ( 6 Tk 2 s + ^ S )'ds 

Jo \ \ [(2/s) 1 q + (Xs) 1 q ] J JO 0 <l<s 

< 8(J 2 V J 6 )A(9-^) A 3+ 4^192(fc 3 ) 2 g 2 A 2 f Je (—L—)ds+ \ E sup |£ / | 2 (6Tfc 2 s + 7 ,)'da, 

Jo V VKr q ) Jo 0 <i<s 

where we have used Lemma 13.31 in the second step. At this point we want to estimate the inverse moments 
of (xt) and to do so we consider the transformation v = x 2 ~ 2q and apply Ito’s formula to get 

(4.3) v t = Vo+[ I (1 - 2q)(l - q)(k 3 ) 2 + 2(1 - q)ki(v s )^ -2(1- q)k 2 v s \ ds+ [ 2k 3 (l - q) y/v^dW s , 

JO " - v -' '- v -' S - v -' JO S - v -' 


K 3 


K 0 K 1 K 2 

for t G [0, T], where vq = (xo) 2 ~ 2q > 0. Denote the drift coefficient of the process (v t ) by a(v t ) and consider 
the function 


(4.4) 


Cl 7 \ \ I TS 1 (^9 — 1)(A + JCo) 29 1 

a(v) := a{v) — A + K 2 v H-- v , 

(fcl)5i=T 


vW 


2-2 q 


where A > 0. Some elementary calculations show that this function attains its minimum at 
and a(v*) = 0 , thus 

a(v) > A — (K 2 + 77 (A)) v. 

Consider the process (Ct(A)) defined through 

(4.5) C*(A) = Co + [ (A - (K 2 + 77 (A)) C s)ds + [ K 3 ^ s dW s , 

Jo Jo 

for t e [0, T] with Co(A) = u 0 . Process (14.51) is a square root diffusion process and when ^4 — 1 > 0 or 

(4.6) A > 2(1 — q) 2 {k 3 ) 2 , 

is a CIR process which remains positive if Co(A) > 0. By a comparison theorem TT] Prop. 5.2.18] it holds 
that v t > Ct(A) > 0 a.s. or (u t ) _1 < (Ct(A )) -1 a.s. or equivalently (x t ) 2q ~ 2 < (Ct(A )) -1 a.s. The inverse 
moment bounds of (Ct(A)) follow by [ 8 ) Rel. 3.1] 

A 


sup E(Ct(A)) p < 00 , for p > -2—2 
te[o,T] 


by choosing big enough A and particularly such that (14.61) holds strictly. Therefore, 

(4.7) E sup Iftj 2 < 8(J 2 V J 6 )A^“ a) A 3 + f E SU p \£ l \ 2 {QTk 2 s + "fs)'ds. 

0 <£<r Jo 0 <l<s 

Relation (14.71) for r = 77 implies 

E sup ( £ t ) 2 < 8(J 2 V J 6 )A (9_ 5) A i + f E sup (£i) 2 ( 6 Tk 2 s + 7 s )'ds 

0 <t<Ti Jo 0 4Z<s 

rl 

< 8(J 2 V J 6 )A (9 -5) A 3 + / e sup ( £ T] ) 2 du 

Jo 0 <j<u 

< 8(J 2 V J 6 )e / A (9 -3 )A 3, 


(4.8) 
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where in the last step we have used Gronwall’s inequality. Using again relation m for t = T and under 
the change of variables u = GTfos + y s we get, 


r0k,2T z -\-^t 


E sup (£ Tj ) 2 du 

0<j<u 


E sup (f t ) 2 < 8(J 2 V J 6 )+ f 
0 <t<T Jo 

< 8(J 2 V J 6 )A (9 "5) a 3 + f e( sup (I {6k2T ^+ 1T > u }£r j ) 2 '] du 

Jo Vo <j<u J 


< 


rbk2l rco / \ 

/ E sup (£ Tj ) 2 du + / P( 6 /c 2 T 2 + 72 -> u)E ( sup (fr,,) 2 |{ 6 fc 2 T 2 + 77 - > u} ) du 

Jo 0 <j<u J&k 2 T 2 \0<j<u J 

+ 8 (J 2 V Je)A^~^ A i 

pOO 

< 8 (J 2 VJ 6 )e 6 i 2 T 2 A («->4 / P( 7 t>u)E sup (f T ,) 2 dit + 8 (J 2 V J 6 )A( ,? -^ a i 

JO 0<j<n 

poo 

< 16(J 2 V J 6 )e 6 /C 2 T 2 A(«-3) A 3+ 8 (J 2 V / P( 7t > u)e“du, 

Jo 

where in the last steps we have used (Id.81) . We proceed by showing that u i-A P( 7 t > u)e u £ £ 1 (M + ). It 
holds that 

(4.9) P( 7t > u) < e _e “E(e C7T ), 

for any e > 0 by Markov inequality. The following bound holds 

fT 192 (fc 3 )V 


7 t = 


'0 [(y s ) 1-9 + (a^) 1 " 9 ] 


! „ r 1 

- 2 ds < 192 (k 3 ) 2 q 2 / ( x s ) 2q 2 ds, 

1_9 1 Jo 


thus 

(4.10) 


E(e £7T ) < E ^ e el92 ( fc 3) 2 9 2 / 0 T U a ) 25 2 *^ ? 


where —1 < 2g — 2<0. It remains to bound the exponential inverse moments of ( Xt ) defined through the 
stochastic integral equation m • Exponential inverse moments for the CIR process are known |141 Th. 3.1] 
and are given by 


(4.11) 


Ee 


«/o*(C.(A))- 1 «J- < C HK (Co)-^ X) -^ X)2+8 ^\ 


for 0 < 5 < — 1^ JJg*- =: ^(A) 2 -^-, where the positive constant Chk is explicitly given in P3J Rel. (10)] 

depends on the parameters fc 2 , k 3 , T, q, but is independent of Co- Thus the other condition that we require 


for parameter A is 
(4.12) 


A > 2(1 - q)V2S(k 3 ) + 2(1 - q) 2 {k 3 ) 2 . 


When (14.121) is satisfied then (14.61) is satisfied too, thus there is actually no restriction on the coefficient 6 in 
(14.111) since we can always choose appropriately a A such that (14.121) holds. Relation (14.101) becomes 

( 413 ) E(e £7T ) < E ^ e el92 ( fe 3)V fj(.v s )~ 1 ds^j < ]g ^ e ei92(fc 3 )V Sj (C»(a)) _1 ^ _ 

We therefore require that 

(4.14) 192(fc 3 )Ve < (^(A)) 2 ^ 


and can always find a e > 1, such the above relation holds by choosing appropriately A as discussed before. 
Relation (14.131) becomes 




and therefore 


E(e^ T ) < C H k((o)- — , 

P(7t > u) < C HK (x o) (1 -« )y(A) e- £ “, 
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where A is chosen such that (14.141) holds with e > 1 . We conclude 

pOO 

E sup (£ t ) 2 < 16(J 2 V J 6 )e 6 /C 2 T 2 A(«-5) A 3+8(J 2 V J 6 )C' ffif (x 0 ) ( 1 "' 3 )y(A) A(«-^ A J / e^~ e)u du 

0 < t<T Jo 

< C- A (9 “5 )A 5, 

by choosing e > 1, where C = C{k\,k 2 , £ 3 , T, q, e) := 8 ( J 2 V J§) ^ 2e 6k2T2 + , is as given 

statement of Theorem 12.21 

5. Alternative approach improving the rate of convergence. 


Proof of Theorem 12.31 We will discuss the proof and highlight the differences since one can follow the proofs 
of Theorems 12.II and 12.21 First of all note that Lemmata 13.11 [3721 and ITTTTT1 still hold, i.e. the moment bounds 
and error bounds of ( yf D ), as well as the moment bounds involving the auxiliary process ( h t ) are true. The 
error £ t := ht — x t now reads as 

(5.1) ht — x t = / (fe{y§,ys) - fe(x 8 ,x 8 )) ds + / {g(y 8 ,y s ) - g(x a ,x 8 )) dW s . 

Jo Jo 

As regards the first part of Theorem 12.31 we now have that 

(5.2) 

0 <t<T 

for any m > 1, where e m = e - m ( m + 1 )/ 2 an d J 3 as stated in Proposition 13.51 The bound (15.21) follows in the 
same lines as in the proof of Proposition 13.51 where now (13.161) becomes 

(5.3) | g{ys,ys) - g(x 8 ,x s )\ 2 < 2(fc 3 ) 2 ( y a \y$ - y s \ 2q ~ l + | h s - y 8 \ + \h 8 - z s |) 
and the term (13.171) disappears. Moreover, we have that 

c t 


/ A 9-5 1 \ 

sup E|£ t | < J 3 -+ 2(fc 3 ) 2 T— e fc2T , 

\ me™ to / 


(5.4) 


E sup \£ t \ < . . 

0 < t<T y /] iL ( A)~ L 


where C e is independent of A and given by C t := 32y ^(fc 3 ) 4 T 2 e 6T fe2+fe2T , where e is such that 0 < e < q— 

The bound (15.41) follows in the same lines as in the proof of Proposition 13.61 where now M t := f*{g(jju, yz ) — 
g(x u ,x u ))dW u and 

E|M t | 2 < 2 (k 3 ) 2 f E (y s \y 8 - y^- 1 ) ds + 2 (fc 3 ) 2 [ - 2 / s |ds + 2 (fc 3 ) 2 / E\£ a \ds, 

Jo Jo Jo 


implying 


E sup \£ t \ 2 < ^(fc 3 ) 2 TJ 3 e fc2T+6T2fe2 ^-^ +32(fc 3 ) 4 T 2 e fe2T+6T2fc2 —, 
o<e<T 3 me m m 


which in turn gives 


1 


E sup \£ t \ 2 < 32(fc 3 ) 4 T 2 e fe2T+6T k2 - _ , 

0 < t<T VlnA"T 

where now we take 0 < e < q — i. 

As regards the second part of Theorem 12.31 we follow Section [2 where now we use the process 

rt 128(fc 3 )V 


7(f) := f 

Jo 


0 [(ye) 1 -* + 


ids 


and the estimate 


\g(ys,Vs) - g{x s ,x s )\ 2 < 2 (k 3 ) 2 y s \y s - y s \ 2q 1 + ^ (| h 8 - y s | 2 + \£ a \ 2 ) (7,)', 
to get the upper bound 

E|M t | 2 < + l -j^ y/E\h a - y s \WmTs)') 2 ds E\£ s \ 2 ( ls )'ds, 
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which in turn implies first 


and then 


E sup (£ t f < ^-J 6 e l A<«-*> 

0 <t<Tl 3 


E sup (f t ) 2 < ^J 6 e 6fc2T2 / P( 7 t > u)e u du. 
o <t<T 3 6 J o 

Finally, in order to bound E ^ e el28 ( fc 3) 2 9 2 / 0 T ( x s) 2q 2 ds ^ , we re q U j re that 


7^2 

2 A 3 


(5.5) 128(fc 3 )Ve < (^(A)) z ^ 

and can always find a e > 1, such that (15.51) holds yielding 

E sup (f t ) 2 <C-A^~^ 


0 <t<T 


where C = C(ki,k 2 ,k 3 ,T,q,e) := (2e 6k2T2 + , is as given in statement of Theorem 

□ 


6. Numerical Experiments. 


We discretize the interval [0,T] with a number of steps in power of 2. The semi-discrete (SD) scheme is 
given by 
( 6 . 1 ) 


v SD - 
Vtn + l ~ 


'yt„ 


i - 


k 2 A 


fci A 


(fc 3 ) 2 A 


1 + k 2 8 A / 1 + k 2 9A 4(1 + k 2 9A) 


;(yt n ) 2q 


-l 


2(1 + k 2 e A) 


{yt n ) q ~ 2 AW n 


for n = 0,..., N — 1, where A W n := Wt n+1 — Wt n are the increments of the Brownian motion which are 
standard normal r.v’s. 

The ALF (Alfonsi) scheme [TJ Sec. 3] is an implicit scheme which requires solving the nonlinear equation 

(6.2) Y n+1 = y tn + (1 - q) (^(K+i)^ - k 2 Y n+l - ^Ml(y n+1 )-^ A + k 3 (l - q)AW n , 

and then computing yfj+i = (Y n+ 1 )~. The estimation of Y n+1 in (|6.2I) can be done for example with 
Newton’s method, but requires a small enough A0 We also consider a scheme recently proposed in [12] using 
again the SD method, but in a different way, 

(6.3) y? n +i(q)= (yt n {l - k 2 A) + fciA - 9 A {y tn ) 2q ~ 1 ^ + k 3 (l-q)AW n 

for n = 0,..., AC — 1. Note the similarity in the expressions of (16.31) and the SD scheme (16.11) proposed here. 
This is not strange, because they both rely in the same way of splitting the drift coefficient. In particular, 
in the explicit HAL scheme, the following process is considered 

yf AL (q) = yt n + 7i{yt n )-A+[ f 2 {y s )ds+ [ sgn(z s )g(y s )dW s , 


for t £ (t n ,t n - )_i] with yo = xq a.s. where now 

(6.4) f(x) = ki- k 2 x - < 77aL x 2 <i- 1 


+ < 77±L x ‘ 2 i 1 . ff ( x ) = k 3 x q 


h( x ) 


and 

(6.5) 


h{x) 
1-9 


Zt = 


(yt n { 1 - k 2 A) + k!A - ^^(y t J 29 " 1 ) + k 3 ( 1 - q)(Wt - W t J. 


*^In the CIR case, i.e. when q = 1/2 116.21) simplifies to a solution of a quadratic equation. 
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A comparison with (|2.2) ) and (|2.3p shows / 2 (:r) = 2qf 2 (x) and g(:r) = for 0 = 0. We write (16.41) again 

as 


(6.6) y? AL {q) = y tn + (4, - fc 2 y t „ - . A + 

and the process (16.61) is well defined when 


gfe ) 2 


(y s ) 2<? 1 ds + k 3 Jt sgn{z s ){y s ) q dW s 


(6.7) 


(fc 3 ) 2 < — fci and A < 


2fc 2 + g(fc 3 ) 2 


The reader can compare again with (12.41) for 9 = 0. Solving for y t , we end up with y^ AL (q) = \z t \ 1 ^ q . The 
main result in [12] is 


E|yf - AL - x t \ 2 <C-A 2q{q ~^, 


when (16.71) holds, implying a rate of convergence at least q(q— ^) which is bigger than the rate of convergence 
of the SD scheme proposed here which is at least \{q — k) ( see Th l2.3D . 

We also consider two more linear-implicit schemes that were stated in the introduction and discussed 
in Appendix [A] Namely, we compare with the balanced implicit method (BIM) with appropriate weight 
functions to guarantee positivity m Th. 5.9]), which reads 

o\ bim/ a _ yt n + fcA + k 3 (y tn ) q (AW n + \AW n \) 

Vtn+iW l + k 2 A + k 3 (y tn )i-i\AW n \ 

and the balanced Milstein method (BMM) with the suggested weight functions [lB] Th. 5.9] that is given by 

ytn + (h + (0 - l)k 2 y tn )A + k 3 (y tn ) q AW n + sH^l(y tn )^(AW n ) 2 


(6.9) 


*Ci (ff) = 


1 + Qk 2 A 


gfe ) 2 


\yt n \ 2q ~ 2 A 


We take the relaxation parameter 0 to be 1/2 as recommended in [lB] Rel. 5.10]. 

We aim to show experimentally the order of convergence for the above positivity preserving methods 
for the estimation of the true solution of the CEV model ED, i.e the semi-discrete methods SD method 
(ED and the HAL scheme (16.31) . as well as the implicit ALF scheme (16.21) and the linear-implicit schemes 
BIM and BMM. The choice of the parameters is the same as in [L5] Fig. 6] with k 3 = 0.4. In particular 
(x 0 ,ki,k 2 ,k 3 ,q,T) = (^, ^, 1,0.4, |, 1). 

Furthermore, we would also like to reveal the dependence of the order of the semi discrete methods on 
q , i.e. we want to verify our theoretical results and in particular the order shown in Theorem 12.31 We take 
the level of implicitness of SD method (16.11) to be 9 = 1, i.e. we consider the fully implicit scheme. We also 
discuss about the fully explicit scheme, i.e. when 0 = 0, but also an intermediate scheme 9 = 1/2, in Section 

El 

We want to estimate the endpoint C 2 — norm e = \/E|y( A l(T) — xt | 2 , of the difference between the nu¬ 
merical scheme evaluated at step size A and the exact solution of ED- For that purpose, we compute M 
batches of L simulation paths, where each batch is estimated by e) = XAi I— y\ r ^\T )\ 2 and the 
Monte Carlo estimator of the error is 


( 6 . 10 ) 


\ 


M L 




(ref). 


ML 


j=l i= 1 


and requires M ■ L Monte Carlo sample paths. The reference solution is evaluated at step size 2 14 of the 
numerical scheme. For the SD case, we have shown in Theorems imio and El that it strongly converges 
to the exact solution. We simulate 100 • 100 = 10000 paths, where the choice for L = 100 is as in I2B p-118]. 
The choice of the number of trajectories M ■ L = 10 4 is also considered in [26] Sec.5] where a fundamental 
mean-square theorem is proved for SDEs with superlinear growing coefficients satisfying a one-side Lipschitz 
condition, but unfortunately it is not positivity preserving. Of course, the number of Monte Carlo paths has 
to be sufficiently large, so as not to significantly hinder the mean square errors. 

We plot in a log 2 — log 2 scale and error bars represent 98%—confidence intervals. The results are shown 
in Table [T]and Figure [Q Table [T] does not present the computed Monte Carlo errors with 98% confidence, 
since they were at least 9 times smaller that the mean-square errors. 
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FIGURE 1. Convergence of fully implicit SD, HAL, ALF, BIM and BMM schemes applied to 
SDE (12.11) with parameters ( xo , fci, k 2 , & 3 , q , T ) = (yg, 1, 0.4, |, 1) with 32 digits of accuracy. 



Step A 

98% SD-Error(0 = 1) 

98% HAL-Error 

98% ALF-Error 

98% BIM-Error 

98% BMM-Error 

2~° 

0.0351607273610989 

0.0356800652730259 

0.0443426636440914 

0.0332754182024868 

0.0358194304596254 

2~' 

0.0350820054172969 

0.035033855267215 

0.0445388652736149 

0.0335414964013289 

0.035736599920736 

2 -y 

0.0345654286067145 

0.0351090011639902 

0.0202607612841043 

0.0329881367906935 

0.0351924591356631 

2 -ii 

0.0332045173200957 

0.0337092293915926 

0.0195315499600022 

0.0317244587795127 

0.0337311936077772 

2 -is 

0.0250782316352445 

0.0249779540239864 

0.0146002661407415 

0.0249983526384181 

0.025291682868944 


TABLE 1. Error and step size of fully implicit SD, HAL,ALF, BIM and BMM scheme for (12.11) 
with (xo, fei, fc 2 , fc 3 , q, T) = (T, A, 1, 0.4, |, 1) and 32 digits of accuracy. 


In Tabic [2] we present the computational times0of fully implicit SD, HAL, ALF, BIM and BMM, for the 
same problem. Figure [2] shows the relation between the error and computer time consumption. As one can 
see from Tabic [2] the CPU times for ALF are at least 1000 times bigger than the other schemes, thus we 
choose in Figure [2] to restrict our attention to the rest of the methods. 


q = 0.75 

Step A 

Implicit SD 

HAL 

ALF 

BIM 

BMM 

Time/Path(in sec) 

2 -b 

0.000013 

0.0000164 

0.0221883 

0.0000174 

0.0000196 

Time/Path(in sec) 

2“' 

0.0000422 

0.0000558 

0.0841705 

0.0000584 

0.0000657 

Time/Path(in sec) 

2 -M 

0.0001586 

0.0002137 

0.2453943 

0.0002207 

0.0002482 

Time/Path(in sec) 

2 -ii 

0.0006243 

0.0008437 

0.9768619 

0.0008703 

0.0009795 

Time/Path(in sec) 

2 ~ Li 

0.0024975 

0.0033977 

3.9096332 

0.0034785 

0.0039143 


TABLE 2. Average computational time (in seconds) for a path, for different discretizations, for all 
considered positivity preserving methods for the mean-reverting CEV process (EH) with q = 0.75. 


Finally, Table [3] presents the exact values of order of convergence for SD, HAL,ALF, BIM and BMM 
produced by linear regression with the method of least squares fit, in the case one considers in one case 3 
points with steps A = 2~ 9 ,2 -11 ,2 -13 and in other case all 5 points including A = 2 -5 , 2 - '. 

We show below, in Tabled the C 2 —distance between our proposed method and the other methods for 
the numerical approximation of m- We work as before and estimate the distance 


( 6 . 11 ) 


d{G,H) = 


\ 


1 M L 
1 (A,G) 


ML 




i=i i =i 

between method G and H , by considering sufficient small A, and in particular for A = 10 —2 , 1CD 3 ,10 -4 . 


4 We simulate with 3.06GHz Intel Pentium, 1.49GB of RAM in Matlab R20146 Software. The random number generator is 
Mersenne Twister. The evaluated times do not include the random number generation time, since all the methods we compare, 
involve the same amount of random numbers. 
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FIGURE 2. Strong convergence error of the mean-reverting CEV process (12.111 as a function 
of CPU time (in sec) using positivity preserving schemes SD,HAL, ALF, BIM and BMM with 
( xo , k i, k 2 , k 3 , q, T) = (-h, 100, 0.05, yg, 1, 0.4, §, 1), and 32 digits of accuracy. 



Time t (in sec.) 


q 

N° of Points 

SD(6> = 1) 

HAL 

ALF 

BIM 

BMM(0 = i) 

0.75 

5 

0.053 

0.054 

0.220 

0.045 

0.054 

3 

0.116 

0.123 

0.118 

0.100 

0.119 


TABLE 3. Order of convergence of SD with 6=1, HAL, ALF, BIM and BMM approximation of 
ED with (x 0 ,ki,k 2 ,k 3 ,q,T) = (A, T, 1,0.4, §,1). 


Step A 

98% - d(SD, HAL) 

98% - d(SD, ALF) 

98% - d(SD, BIM) 

98% - d(SD, BMM) 

Icr 2 

0.000572742828312345 

0.0716139598501867 

0.0038372799928164 

0.000531186647260291 

10 -d 

0.000157661207152358 

0.0286630459661306 

0.00134601166557609 

0.000156449754950229 

icr 1 

0.000049813956973071 

0.0283116960890337 

0.00044483884745991 

0.000049771403122846 


TABLE 4. The C 2 —distance between all the considered numerical schemes applied to SDE ED 


with parameter set (zo, fci, k 2 , k 3 , q, T) = (^, i ; 0.4, |, 1). 


Finally, we examine the behavior of SD w.r.t the parameter q. We want to examine the impact of q in 
the order of convergence and verify our theoretical results and in particular Theorem 12.31 Table [5] shows the 
order of SD w.r.t q. 


q 

Order of Fully Implicit SD 

0.6 

0.110 (0.052) 

0.7 

0.118 (0.0525) 

0.8 

0.116 (0.05) 

0.9 

0.121 (0.053) 


TABLE 5. Order of convergence of SD(# = 1) approximation of (12.11) with (xo, ki, k 2 , k 3 , T) = 
(T, i ) 0.4,1) for different values of q. 


The following points of discussion are worth mentioning. 

• The performance of all methods, as shown in Table [l] and Figure [lj implies, in terms of error 
estimates, that the implicit ALF scheme performs better, for values of discretization steps A > 2~ 9 . 
All the other methods, i.e. the semi discrete SD and HAL, and the BIM and BMM have a similar 
behavior for all A’s in the sense above as Figure [U shows. The similarity of SD, HAL, BIM and 
BMM is also indicated in Table 2) where we see how close they are w.r.t. the C 2 — norm, and in 
Table [3] where the convergence order is considered. Nevertheless, Table ID also shows that in order to 
get an accuracy to at least 2 decimal digits, which in practice may be adequate concerning that we 
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want for example to evaluate an option and thus our results are in euros, there is no actual harm in 
choosing whatever of the above available methods. We may then choose the fastest one, as will be 
discussed later on. 

• We see that the strong order of convergence of implicit SD for problem m is at least l/2(q — 1/2) = 
1/8, as shown theoretically and presented in Tabic [3] We also see that all methods converge with 
similar orders and the theoretically rate 1 of the ALF method [I] does not hold for these A’s. Thus, 
again we see that the rate in practical situations does not necessarily matter, if one has to compute 
with very small A’s to achieve it. Moreover, we present in Table [6] the performance of the explicit 
SD method and see that it is very close to the implicit, which is of course natural to happen. 


Step A 

98%-SD-Error(6> = 0) 

Order 

2~° 

0.034005669022654 

0.113 

(0.047) 

2~‘ 

0.0344244107574687 

2 -a 

0.0342415044537196 

2 -ii 

0.0331273071237492 

2 T7T3 

0.0250195459354763 


TABLE 6. The performance of fully explicit SD scheme (16.111 applied to SDE (12.11) with parameter 
set (x 0 ,k 1 ,k 2 ,k 3 ,q,T) = (A, A, 1,0.4, |, 1) 


• In practice, the computer time consumed to provide a desired level of accuracy, is of great importance. 
Especially, in financial applications, a scheme is considered better when except of its accuracy, it 
is implemented faster. As mentioned before, the SD method as well as the HAL method performs 
well in that aspect, compared to the implicit ALF method, which requires the estimation of a root 
of a nonlinear equation in each step and is therefore time consuming. This is presented in Table [2] 
and Figure [2] which illustrates the advantage of the semi-discrete method SD, performing slightly 
better than HAL and BMM, better than BIM, and of course a lot better compared with ALF (over 
1000 times quicker to achieve an accuracy of almost 2 decimal digits.) Moreover, the explicit SD, 
performs slightly better in that aspect, as shown in Table [7] 


q = 0.75 

Step A 

Fully Explicit SD (Implicit) 

Time/Path(in sec) 

2 ~ b 

0.000013 (0.000013) 

Time/Path(in sec) 

2 ~‘ 

0.0000411 (0.0000422) 

Time/Path(in sec) 

2 -y 

0.0001545 (0.0001586) 

Time/Path(in sec) 

2" 11 

0.0006048 (0.0006243) 

Time/Path(in sec) 

2 -13 

0.0024319 (0.0024975) 


TABLE 7. Average computational time for a path (in seconds) for fully explicit SD method for q = 0.75. 


• A negative step of a numerical method appears when the computer-generated random variable ex¬ 
ceeds a certain threshold, which tends to increase as the step size A decreases. Thus, the undesirable 
effect of negative values that are produced by some numerical schemes (such as the explicit Euler 
(EM) and standard Milsten (M) ), tends to disappear, since after a certain small step size, the thresh¬ 
old exceeds the maximum standard normal random number attainable by the computer system. 


7. Approximation of stochastic model (E3D- 

So far we have focused on the process (V t ), which is one part of the two-dimensional system (11.11) . 
Nevertheless, it can be treated independently, since the only way that it interacts with the process ( S t ) is 
through the correlation p of the Wiener processes. First we apply Ito’s formula on In (St) to get, 

In S t = In S 0 + [ pdu j ( V u ) 2p du + [ (V u ) p dW u , t G [0, T\. 

Jo ^ Jo Jo 


(7.1) 
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Then, we consider two different schemes for the integration of (0)0 The first is the EM scheme which 
reads 

(7-2) InS^" = lnS tn + /xA - i(U t J 2p A + (V tn ) p AW n , 

has strong convergence order 1/2 and is easy to implement. The second scheme, which is based on an 
interpolation of the drift term and an interpolation of the diffusion term, considering decorrelation of the 
diffusion term, including a higher order Milstein term m Sec.4.2], is denoted IJK and is given by [15, 
Rel.(137)] 


In S"* = In S tn + pA - ± ((V tn ) 2p + (V tn+1 ) 2p ) A + p(V tn ) P AW n 

(7-3) ((V u y + ( V tn+1 Y ) (A W n - P AW n ) + \ppk 3 (V t J’+P" 1 ((A W n ) 2 - A) . 

We therefore consider the EM scheme (17.21) combined with SD (16.111 . the IJK scheme (17.31) combined with 
SD (16.11) and compare with the case where the stochastic variance (p = f) is integrated with BMM scheme 
(16.91) . for three different correlation parameters, p = 0,p = —0.4 and p = —0.8 with So = 100, p = 0.05, as 
in EE3 Sec.5]. We present in Tables [51 l9l and ITOl and Figures [3J [4] and [3 the errors, in the sense of distance 
(16.111) . for all the above considered ways of numerical integration of process (St), for different step sizes, as 
well as the average computational time (in seconds) consumed for each discretization. 


Step A 

EM&;SD-Error(0 = 0.5) 

IJK&SD-Error(<9 = 0.5) 

EM&BMM-Error(© = 0.5) 

IJK& BMM-Error(© = 0.5) 

2~° 

26.901 (0.0000261) 

26.901 (0.0000159) 

26.891 (0.00002) 

26.890 (0.0000294) 

2~‘ 

27.288(0.0000919) 

27.288(0.0000492) 

27.277 (0.0000676) 

27.277 (0.0001043) 

2 -y 

27.298(0.0003595) 

27.297(0.0001843) 

27.289 (0.0002610) 

27.288 (0.0004081) 

2 -11 

25.057(0.0014255) 

25.058(0.0007309) 

25.051 (0.0010309) 

25.051 (0.0016191) 

2~ 13 

19.441 (0.0057322) 

19.441 (0.0028928) 

19.442 (0.0041177) 

19.442 (0.0064721) 


TABLE 8. 98%—Error, step size and average computational time of numerical integration 
of process (St) using log-Euler or IJK method with SD or BMM scheme for (11.11) with 
(*o, So, p, ki, fe, & 3 , q, T) = (A, 100, 0.05, A, 1,0.4, |, 1), correlation p = 0 and 32 digits of accuracy. 


FIGURE 3. Strong convergence error of the financial underlying process (St), as a function 
of CPU time (in sec) using log-Euler or IJK method with SD or BMM scheme for CD with 
(xo, So, p, k\, & 2 , & 3 , q, T) = (A, 100, 0.05, A, 1, 0.4, |, 1), correlation p = 0 and 32 digits of accuracy. 



Figures [3l 0] and [5] indicate that in all cases the favorable choice is to integrate (S t ) using IJK method 
combined with the SD scheme for (V*) in model (11.11) . The above combination seems to be the better one, 
w.r.t. CPU time, for every correlation coefficient considered. 

U he reason for not considering other schemes such as the two-dimensional Milstein is that they generally are time consuming, 
since they involve additional random number generation for the approximation of double Wiener integrals. 
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Step A 

EM&SD-Error(0 = 0.5) 

IJK&SD-Error(0 = 0.5) 

EM&BMM-Error(© = 0.5) 

IJK& BMM-Error(0 = 0.5) 

2~° 

26.382 (0.0000266) 

26.331 (0.0000161) 

26.372 (0.0000202) 

26.324 (0.00003) 

2~‘ 

26.448(0.0000951) 

26.396(0.000005) 

26.439 (0.0000691) 

26.389 (0.0001081) 

2 -y 

25.951 (0.0003631) 

25.909(0.000184) 

25.944 (0.0002606) 

25.904 (0.0004131) 

2 -ii 

24.540(0.0014506) 

24.494(0.0007355) 

24.531 (0.0010378) 

24.486 (0.0016495) 

2 —ld 

18.738(0.0060748) 

18.749(0.0030185) 

18.735 (0.0042868) 

18.747 (0.0068395) 


TABLE 9. 98%—Error, step size and average computational time of numerical integration 
of process (St) using log-Euler or IJK method with SD or BMM scheme for (11.11) with 
(xo, So, p, fci, & 2 , & 3 , q, T ) = (-b, 100, 0.05, jb, 1, 0.4, |, 1), correlation p = —0.4 and 32 digits of ac¬ 


curacy. 


FIGURE 4. Strong convergence error of the financial underlying process (St), as a function 
of CPU time (in sec) using log-Euler or IJK method with SD or BMM scheme for E3 with 
(xo, So, Ah ki, & 2 , &3, q, T) = (-^, 100, 0.05, 1,0.4, f, 1), correlation p = —0.4 and 32 digits of 

accuracy. 




\\ v. 
\ 


q x Ko 
\ .. 


- * - log-Euler & SD 

- e - IJK & SD 

X log-Euler & BMM 
O - IJK& BMM 


3 4 

Time t (in sec.) 


Step A 

EM&;SD-Error(0 = 0.5) 

IJK&SD-Error(<9 = 0.5) 

EM&BMM-Error(© = 0.5) 

IJK& BMM-Error(© = 0.5) 

2~° 

25.552(0.0000263) 

25.455(0.0000159) 

25.541 (0.0000199) 

25.449 (0.0000296) 

2~' 

25.670(0.0000932) 

25.569(0.0000494) 

25.659 (0.0000678) 

25.564 (0.0001059) 

2 -y 

25.217(0.0003622) 

25.137(0.0001835) 

25.208 (0.0002595) 

25.132 (0.0004111) 

2 -n 

23.743(0.0014407) 

23.711 (0.0007306) 

23.734 (0.0010307) 

23.707(0.0016376) 

2 -I3 

18.082(0.005871) 

18.316(0.0029312) 

18.078 (0.0041637) 

18.312 (0.0066239) 


TABLE 10. 98%—Error, step size and average computational time of numerical integration 
of process (St) using log-Euler or IJK method with SD or BMM scheme for (11.11) with 
(xo, So, p,k\,k 2 ,ks,q,T) = (T ) IQO, 0.05, 1, 0.4, §, 1), correlation p = —0.8 and 32 digits of ac¬ 


curacy. 


8. Conclusion. 

In this paper, we exploit further the semi-discrete method (SD), originally appeared in [TO] , to numerically 
approximate stochastic processes that appear in financial mathematics and are meant to be non negative. 
In [13] we examined the Heston 3/2—model, that is a mean reverting process with super-linear diffusion, 
described by a SDE of the form m with q = 3/2. Now, we deal with SDEs with sub-linear diffusion 
coefficients of the type (xt) q with 1/2 < q < 1. These kind of SDEs, called mean reverting CEV processes, 
appear in stochastic models, where they represent the instantaneous volatility-variance of an underlying 
financially observable variable. We prove theoretically the strong convergence of our proposed SD scheme, 
revealing the order of convergence. The resulting polynomial rate in Theorem l2.3l may not appear appealing 
at first sight, because of its low magnitude. Nevertheless, as it is shown in the numerical experiment section, 
where a comparative study is presented between various positivity preserving schemes, the SD method seems 
to be the best w.r.t. CPU time consumption. The advantage of the SD method here is that although implicit, 
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FIGURE 5. Strong convergence error of the financial underlying process {St), as a function 
of CPU time (in sec) using log-Euler or IJK method with SD or BMM scheme for dJ with 
(xo, So, p, ki, k 2 , k 3 , q, T) = (jb, 100, 0.05, yg, 1, 0.4, f, 1), correlation p = —0.8 and 32 digits of 
accuracy. 



has an explicit formula and thus requires fewer arithmetic operations and consequently less computational 
time. Moreover, our method can cover cases where GU has time varying coefficients, i.e. k\{t), k^t), k^(t). 

We also treat the whole two-dimensional stochastic volatility model m- In order to do that, we actually 
integrate the process ln(St) which satisfies a SDE of the form (17.111 and in the end transform back for (St). 
We only consider two different schemes for the integration of ln(5*), namely the Euler Maruyama (EM) 
scheme, which is easy to implement and the IJK scheme [15] Rel.(137)] which is shown to be the most 
efficient method, robust and simple as EM [T5] - We do not apply other two-dimensional schemes, such as 
for example the Milstein scheme, since they are in general time consuming, as they involve approximations 
of double Wiener integrals which require additional random number generation. We therefore combine the 
EM scheme with SD (o & mi the IJK scheme with SD ( (17.31) & (16.11) ) and compare with the case 
where the stochastic variance (p = J) is integrated with BMM scheme (16.91) . for three different correlation 
parameters, p = 0,p = —0.4 and p = —0.8 with So = 100,/x = 0.05, as in [T5] Sec.5]. The combination IJK 
with SD seems to be the most favorable w.r.t. CPU time, for all the cases. 
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Appendix A. Some numerical schemes for the integration of the variance-volatility 

PROCESS (Vj). 

We consider a partition of the time interval [0, T] with 0 = to < D < ••• < ijv = T and discretization steps 
A„ := t n+ 1 — t n for n = 0,..., N — 1. Moreover, we denote by A W n := W tn+1 — W tn the increments of the 
Brownian motion. We show in the following subsections some numerical schemes for the approximation of 

(A.l) Vt = V 0 + [ (fci - k 2 v s )ds + f k 3 (V s ydW s , t G [0, T] 

Jo Jo 

and make some brief comments on them. We also denote V n := Vt„ ■ 

Standard Euler-Maruyama scheme. The Euler method, applied to the SDE setting, already appeared 
in the 50's through Maruyama [22] and thereafter there has been an extensive study on numerical approxi¬ 
mations of solutions of SDEs (we just mention [lB] for a recent review on numerical methods for SDEs with 
applications in finance and references therein). 

The explicit Euler-Maruyama (EM) scheme for the process (V)) is given by 

(A.2) V n E ff =V n + (fcr - k 2 V n ) A„ + fc 3 (V„)«AW n , 

for n = 0,..., N — 1. Clearly P(U+i < 0|U > 0) > 0, thus the EM scheme can produce negative values 
with positive probability, or in the notion of [25] we say that (IA.2I) has a finite life time. 
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Standard Milstein scheme. The standard one dimensional Milstein (M) scheme contains some extra 
terms derived by Ito-Taylor expansion m Sec.5], and applied to (V t ) reads 

(A.3) K+i =V n + (fci - k 2 V n )A n + k 3 (V n ) q AW n + \{k 3 ) 2 g(K) 29 " 1 ((A W n ) 2 - A„) , 

for n = 0,..., N — 1 where we have retained terms of order (A„). Again (M) scheme has a finite life time. 


Balanced Implicit Method. The balanced implicit method (BIM) [23] Rel (3.2)] was the first attempt to 
treat the problem of invariance-preserving of specific domains of the underlying process and reads 

(A.4) = V n + (fci - k 2 Vn)An + k 3 (V n ) q AW n + (c°(V„)A n + c\V n )\AW n \) (V n - V n+1 ), 

for n = 0,..., N — 1 where c° and c 1 are appropriate weight functions. The choice c°(a;) = k 2 and c 1 (x) = 
k 3 x q ~ x preserves positivity [16] Sec. 5]. Rearranging the above equation, we get the expression 

Vn + fclA n + k 3 (V n ) q {AWn + \AW n \) 


(A.5) 


V 


BIM 

n+1 


l + fc 2 A n + kaiVn^-^AWn 


Balanced Milstein Method. The balanced Milstein method (BMM), was proposed in [16], for an im¬ 
provement of the BIM in the stability behavior but also both in the rate of convergence. It is given by the 
following linear implicit relation 


t /BMM 
v n +1 


V n + {h - k 2 V n ) A n + k 3 {V n ) q AW n + ^ (fc 3 f q(V n ) 2q ~ 1 ((A W n f - A n ) 
+ {d°(V n )A n + d\V n )((AW n ) 2 - A n )) (Vn - V n+1 ), 


for n = 0, ...,N — 1 where d° and d} are appropriate weight functions. The choice d°(x) = Qk 2 + 
\{k 3 ) 2 q\x \ 2q ~ 2 , where 0 £ [0,1] and d 1 (x) = 0 implies an eternal life time for the scheme [16] Th. 5.9], in the 
sense that P(T4+i > 0|VA > 0) = 1. The step sizes A ra have to be such that A n < y ■ The relaxation 

parameter resembles to the implicitness parameter (0 in our notation). For 0 = 1 there is no restriction in 
the step size, but it is recommended when possible na Rem. 5.10] to take 0 = 1/2. Rearranging with the 
above specifications leads to 


(A.6) 


t rBMM 

v n +1 


V n + (k 1 - (1 - 0)k 2 V n )A n + k 3 {V n ) q AW n + \(k 3 ) 2 qiVn) 2 ^ 1 {AW n ) 2 
l + Qk 2 A n + l(k 3 ) 2 q\V n \ 2 <i- 2 An 


Finally, the proposed semi-discrete (SD) scheme reads 
(A.7) 


V SD — 

v n +1 — 




k 2 A 


1 + k 2 0 A 


k\A 

1 + k 2 6A 


(k 3 ) 2 A 
4(1 + k 2 9A) 2 


(Vn ) 2 ^ 1 + 


2(1 + k 2 6A) 


(V n ) q -*A W n 


Increasing the time horizon T results in an increase of the percentage of negative paths of EM and M. On 
the other hand BIM, BMM and of course SD are not affected by that, since they preserve their positivity on 
any interval [0,T], 
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